3.8 This works on lines and polygons as well!:

We can apply this same logic and format to query and spatially display polygon based layer. For example, we will query the TMDL Watersheds layer (Layer ID 72), searching for all Benthic Impairment TMDLs. We can also convert this data and read it into a data frame for further interpretation:

url <- parse_url("https://gis.deq.virginia.gov/arcgis/rest/services")
url$path <- paste(url$path, "staff/DEQInternalDataViewer/MapServer/72/query", sep = "/")
url$query <- list(where = "IMP_NAME = 'Benthic-Macroinvertebrate Bioassessments (Streams)'",
                  outFields = "*",
                  returnGeometry = "true",
                  f = "geojson")
request <- build_url(url)

Benthic_TMDLs <- st_read(request)


pal <- colorFactor(
      palette = "Set2",
      domain = unique(Benthic_TMDLs$WSHD_ID))

CreateWebMap(maps = c("Topo","Imagery","Hydrography"), collapsed = TRUE) %>%
  addPolygons(data = Benthic_TMDLs, 
              color = ~pal(WSHD_ID),
              fillColor = ~pal(WSHD_ID), 
              fillOpacity = 0.6,
              stroke=5,
              label = ~WSHD_ID )